// Copyright (c) 2005  CRIMERE (France).
// All rights reserved.
//
// This file is part of OPSC; you can redistribute it and/or
// modify it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; version 2.1 of the License.
// See the file LICENSE.LGPL distributed with CG-LAB.
//
// Licensees holding a valid commercial license may use this file in
// accordance with the commercial license agreement provided with the software.
//
// This file is provided AS IS with NO WARRANTY OF ANY KIND, INCLUDING THE
// WARRANTY OF DESIGN, MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE.
//
// Author(s)     : Mercier Jean-Marc
//=========================================================================
#include <set>


#include <src/algorithm/random/BoxMuller.h>
#include <src/algorithm/random/UniformEquallySpacedGenerator.h>
#include <src/algorithm/random/NormalGenerator.h>
#include <src/math/Operators/delta_particle.h>
#include <src/algorithm/EvolutionAlgorithm.h>
#include <src/algorithm/Evolution_Algorithm/CrankNicolson.h>
#include <src/algorithm/Evolution_Algorithm/ForwardEuler.h>
#include <src/math/CauchyProblems/CauchyProblem_traits.h>

/*
	This class initialize the Cauchy problem for multidimensional heat equation.
	The algorithm used is Crank Nicolson.
	The numerical method is particle one.
	The initial data corresponds to a normal distribution.
*/

#if !defined(_CauchyProblem_HeatEquation_)
#define _CauchyProblem_HeatEquation_

#define NUMBER_OF_PARTICLES 9
#define DIM_D 2
#define DIM_P 1

		template <class weight_set, class Time = OS_double>
		class CauchyProblem_HeatEquation_traits : public CauchyProblem_traits< OS::CrankNicolson < Delta_p_t<weight_set> > >
		{
		public :
			CauchyProblem_HeatEquation_traits() 
				: 
				dim_d(DIM_D),
				dim_p(DIM_P)
			{
			};
			virtual void initial_conditions(int nb_part)  = 0;
		protected :

			int dim_d, dim_p;
		};


		template <class weight_set, class Time = OS_double>
		class Cauchy_Delta_p : public CauchyProblem_HeatEquation_traits<weight_set,Time>
		{

			public :
			Cauchy_Delta_p(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};
			virtual void set_dim(int dim_d_)
			{
				assert(dim_d_ > 1);
				dim_d = dim_d_;
				initial_conditions(NUMBER_OF_PARTICLES);
			}
			void initial_conditions(int nb_part)
			{
				particles_list temp(nb_part,dim_d);
				boost::random::BoxMuller<particles_list::particle_type> generator(dim_d);
				generator.fill(temp);
				set_initial_conditions(temp);
				insert_time(0.);
				insert_time(1e+10);
			}

		};

typedef weight_set<particle, OS_double >					OS1d_edges_set_type;
typedef Cauchy_Delta_p<OS1d_edges_set_type>		OS1d_time_grid_type;

		template <>
		class Cauchy_Delta_p< OS1d_edges_set_type,OS_double> : public CauchyProblem_HeatEquation_traits<OS1d_edges_set_type,OS_double>
		{
		public :
			Cauchy_Delta_p(){
				initial_conditions(NUMBER_OF_PARTICLES);
			};

			void initial_conditions(int nb_part)
			{
				dim_p = nb_part;
				particles_list temp(nb_part);
				boost::random::UniformEquallySpacedGenerator local(dim_p);
				boost::random::NormalGenerator<boost::random::UniformEquallySpacedGenerator>	
					generator(0.0, M_SQRT2 *std::sqrt( 0.1 ),local);
				generator.fill(temp);
				set_initial_conditions(temp);

				insert_time(0);
				insert_time(1e+10);
			}
		};

#endif